Orbital magnetic susceptibility of finite-sized graphene 
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We study the orbital magnetism of graphene ribbon in the effective-mass approximation, to figure 
out the finite-size effect on the singular susceptibility known in the bulk limit. We find that the 
susceptibility at T = oscillates between diamagnetism and paramagnetism as a function of ef, 
in accordance with the subband structure formed by quantum confinement. In increasing T, the 
oscillation rapidly disappears once the thermal broadening energy exceeds the subband spacing, 
and the susceptibility x( £ f) approaches the bulk limit i.e., a thermally broadened diamagnetic peak 
centered at ef — 0. The electric current supporting the diamagnetism is found to flow near the 
edge with a depth ~ hv / {2-KksT), with v being the band velocity, while at T = the current 
distribution spreads entirely in the sample reflecting the absence of the characteristic wavelength 
in graphene. The result is applied to estimate the three-dimensional random-stacked multilayer 
graphene, where we show that the external magnetic field is significantly screened inside the sample 
in low temperatures, in a much stronger manner than in graphite. 
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I. INTRODUCTION 

Graphite is known as one of the strongest diamag- 
netic materials among natural substances. 1-3 This prop- 
erty is due to the large orbital diamagnetism related to 
the small effective mass in the band structure, i.e., nar- 
row energy gap between conduction and valence bands. 
The diamagnetic effect becomes even greater in graphene 
monolayer 4-6 which is truly a zero-gap system. 1:7-9 The 
magnetic susceptibility of graphene at zero temperature 
contains a singularity expressed as a delta function in 
Fermi energy ef, which diverges at Dirac point (ef = 0) 
where the two bands stick, and vanishes otherwise. l ' 10-17 
The orbital diamagnetism has been studied for other 
graphene-related materials, such as graphite intercala- 
tion compounds, 18,20-22 carbon nanotube, 23-25 few-layer 
graphenes, 19,26 ' 27 and an organic material having similar 
gapless spectrum. 28 

In this paper, we investigate the orbital diamagnetism 
of a graphene strip with finite width. 29-41 The purpose 
of this work is two-fold: (i) To understand how the delta- 
function singularity of the bulk limit is relaxed in a real- 
istic finite-sized graphene system. In the literatures, the 
orbital susceptibility of graphene nanoribbons was cal- 
culated for the Fermi energies near the Dirac point, 31 ' 44 
while the behavior off Dirac point and the relation to the 
bulk susceptibility is not well understood, (ii) To study 
the diamagnetic current flow on graphene. In the con- 
ventional diamagnetism of metal, we usually expect that 
the current circulates near the surface with a depth of 
the order of the Fermi wave length Xp. In graphene, the 
only characteristic length scale Xp intrinsically diverges, 
and we expect a peculiar current distribution different 
from the conventional system. 

To address above problems, here we calculate the 
orbital susceptibility and the current distribution of 
graphene ribbon with an arbitrary width, in various 
Fermi energies e p and temperatures T, using the effective 
mass approximation. We find that the susceptibility at 
T = oscillates between diamagnetic and paramagnetic 



values in increasing ef, in accordance with the detailed 
subband structure. In increasing temperature, the os- 
cillation rapidly disappears once the thermal broadening 
energy exceeds the subband spacing, and the suscepti- 
bility approaches bulk limit, i.e., a thermally-broadened 
diamagnetic peak centered at Ep — 0, independently of 
the atomic configuration at the edge. We also apply a 
similar analysis to the carbon nanotube, and find a sim- 
ilar oscillation in the susceptibility. 

The electric current supporting the diamagnetism 
spreads entirely in the sample at T = 0, reflecting the 
absence of the characteristic wavelength. In increasing 
temperature, however, the current density tends to lo- 
calize near the boundary with a depth ~ hv/ (27rfcsT), 
forming an edge current circulation. 

The analysis of the spatial distribution of the diamag- 
netic current is useful in studying a graphene stack where 
the diamagnetic current of one layer influences the elec- 
tron motion in other layers. If we take a randomly- 
stacked graphene multilayer, in which the interlayer cou- 
pling is expected to be small, 45-51 the self-consistent cal- 
culation shows that the diamagnetism is much stronger 
than in graphite, and the external magnetic field is signif- 
icantly screened inside the sample in low temperatures. 

The paper is organized as follows. In Sec. II, we briefly 
review effective mass description of electrons in graphene 
ribbon, and formulate the orbital magnetic susceptibility. 
We present the numerical results and detailed discussion 
in Sec. Ill, as well as a similar analysis for the carbon 
nanotube in Sec. IV. We argue the diamagnetism of 
randomly stacked graphene multilayers in Sec. V. The 
conclusion is given in Sec. VI. 



II. FORMULATIONS 

A. Effective mass approximation 

Graphene is composed of a honeycomb network of car- 
bon atoms, where a unit cell contains a pair of sublattices, 
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wave amplitudes are written as 53 

MRa) = e iK - R *F*{R K ) + ^e iK '- R ^Ff(R A ), 

MRb) = -^e lKR *F%{R B )+J K ' R *Fg{R B ). 

(2) 

in terms of the slowly-varying envelope functions 
F A ,F B ,F A ,and F^ . The envelope functions satisfy 
the Schrodinger equation, 1,7-9:52,53 



H Q F(r) =eF{r), 



(3) 
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FIG. 1: Atomic structures of graphene ribbons with (a) zigzag 
boundary and (b) armchair boundary, respectively. Dashed 
circles indicate missing sites beyond the boundary. 



denoted by A and B. Fig.l (a) and (b) show the atomic 
structure of zigzag and armchair graphene ribbons, re- 
spectively, where a and b are primitive translation vec- 
tors of infinite graphene. The lattice constant is given by 
a = |a| « 0.246 nm. For both cases we set y-axis along 
the ribbon, and set x = and L x to the line of missing 
sites nearest from the edge. We define ij as the angle 
between x axis and a, which is 7r/6 for zigzag, and for 
armchair boundary. 

In a tight-binding model, the wave function of 
graphene electron is written as 

Hr) = MRA)<t>(r -R A )+Y, MRB)Hr - Rb), 

Ra Rb 



(i) 



where R A and Rb are the positions of A-sites and B- 
sites, respectively, and <fi(r) denotes the wave function of 
the p z orbital of a carbon atom. 

For states in the vicinity of the Fermi level e = 0, the 
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Ff (r) 



(5) 



where k — — iV and v is the band velocity. 

The electronic states of the graphene ribbon can be 
correctly described by setting the appropriate boundary 
condition to the effective mass Hamiltonian. 34 Now the 
eigenstates are labeled by k y since the system is transla- 
tionally symmetric along y. A wave function of k y and 
the energy e is generally written as 



F(r) 



— e ik v y 



(6) 



where k 2 x = e 2 /h 2 v 2 ~k 2 y , e* e = (k x +ik y ) / \J ' k 2 x + fc 2 , s = 

s/\e\, and A,B,C and D are numbers to be determined 
by satisfying the boundary condition, as we will argue in 
the following. 



B. Zigzag boundary 

In the zigzag ribbon, the boundary condition is given 
by iPa(Ra) — at x — 0, and ^b(-Rb) = at x = L x . 
By using Eq. (2), this is translated to the condition for 
the envelope function as 



*a(0,») = 0, 
F|(L x ,j/) = 0, 

Ff' (0,y) = 0, 
F B i '(L x ,y) = 0, 



(7) 
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FIG. 2: Plots to find solutions of Eq.(9) and Eq.(ll). 

which keeps the states at K and those at K' independent. 
For an eigenstate for the K point, we apply the first two 
lines of Eq.(7) to Eq.(6), to obtain 

A + B = 0, 

s (^A e i ( k * L *+ e ) — g e -'(*« I '«+ 8 )j = o (8) 

To have a solution other than A = B = 0, we require 34 

k. r 



t an 



(9) 



For given ky, we define k n (n — 0, 1, 2, • • • ) as solution 
of Eq.(9) in k x satisfying nir < k n L x < (n + l)ir. Corre- 
sponding eigenstates and the energy are obtained as 



F snky (r) = A T 



A n = 1 - 



ik y y 



2k n L x 



i sin k„x 
s{-l) n+1 sin k n {x-L x ) 



-1/2 



SsUky = ShV^kl + k 2 y, 



(10) 



for n = 0, 1, 2, • • • . For the normalization of the wave- 
function, we assumed the periodic boundary condition in 
y-direction with a large enough period L y . 

Solution k n is obtained by searching for crossing points 
of tank x L x and k x /k y as illustrated in Fig. 2. When 
k y L x > 1, the first solution fco becomes a pure imaginary 
number ikq which satisfies 



tanh kqL x 



(11) 



The wavefunction and the energy then becomes 



F sQky (r) = A! 



A Q = -1 



„ik y y 

y/L X Ly 



/ i sinh KqX 
-s sinh kq(x — L x ) 



sinh2Koia;^ 
2kqL x 
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£ S Qk y = sKvJ-kI + kh 



(12) 



This actually describes the edge state localized at the 
boundary x — and L x giving a nearly flat energy 
band. 29-31 

The eigenenergy s sn k represents the n-th branch of 
conduction (s = +) and valence (s = — ) bands re- 
spectively. Energy band structure of K as a function 
of k y is shown as solid curves in Fig. 3 (a). Eigenstates for 
K' point are obtained similarly, where the energy band 
structure is equivalent to Fig. 3(a) with k y inverted to 
—k y . The flat band of edge states of K and K' are con- 
nected in a wave number away from K or K'. 29-31 



C. Armchair boundary 

In the armchair ribbon, the boundary condition im- 
poses both of iPa{Ra) — and ipB (-Rb) = at each of 
x = and x = L. x . The corresponding conditions for the 
envelope functions are written as 

Fg{0,y)+F$' (0,y)=0, 
Fg(0,y)- Fg'(0,y)=0, 
Ff{L x ,y) + u- 2N Ff\L x ,y) = 0, 
Fg(L Xl y)-u- 2N Fg'(L x ,y) = 0, (13) 

where N — L x /a is the number of honeycomb lattices 
between x = and L x , which can be integer or half- 
integer depending on the position of the edge. 
Applying above conditions to Eq.(6), we obtain 
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(14) 



where a 



-2N 



and A = k x L x . The determinant of ma- 



trix in Eq.(14) should vanish to have a non-zero solution. 
This condition is reduced to 
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7T 

T x 



0,±1,±2,- 



v is an integer (0, ±1) defined by 
2N = 3m + ; 



(15) 



(16) 
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with integer m. The eigenstate and energy are obtained 

as 



Fsnk y (r) 



ik y y 



( 



2\/L x Ly 

£snk y = sHv^Jkl + fc2. 



se 



\se 



i(k n x+8) 
-i(k„x-8) 



(17) 



When v = 0, the energy bands of n — and s = ± 
stick together and thus the system is metallic, while oth- 
erwise a gap opens at zero energy and the system be- 
comes a semiconductor. Energy bands for metallic arm- 
chair ribbon (v = 0) and semiconducting armchair ribbon 
(u = ±1) are shown in Fig. 3(b) and (c), respectively. In 
(c), the labeling n is for the case of v = +1, while n 
becomes —n in v = — 1. 



D. Orbital susceptibility 

To calculate the orbital diamagnctism, we consider a 
graphene ribbon under a uniform magnetic field B per- 
pendicular to graphene plane. We take the Landau gauge 
and set the vector potential as 



A(r) 



0. ^!l-y 



(18) 



The Hamiltonian in presence of the magnetic field is ob- 
tained by replacing k by k + eA/(hc), as 



H = U a + 5H, 5% = -v y A y , 



(19) 



where c is the light velocity, and 



1 dU 



/0 -i 0\ 

i 

i 

\0 -i Oj 



(20) 



case in armchair ribbons. This is obvious from at the 
matrix clement of j y between two eigen states F and F' , 

(F'\j y (r)\F) = ie v [F> A K (ryF«(r) F' B K (vf F% (r)] . 

(23) 

In the wavefunction of zigzag ribbon, Eqs. (10) and (12), 
the component F% is zero at x = 0, and Fg is at L x , so 
that Eq. (23) vanishes at the both edges. 

The current density on xy-plane is related to the local 
magnetic moment m(r) in z-direction by 



Jx 



dm 
dy ' 



dm 
dx ' 



(24) 



In the present case, m(r) depends only on x so that the 
total magnetization per area is 



M 



1 



L X L V j a 



1 



m(x)dxdy 



cL x Jo V ^ 
The magnetic susceptibility is written as 



j y (x)dx. 



(25) 



r M 
X = lim — = 

b^o B cL x L y 



5>-> E » 



(26) 



We calculate Eqs. (22) and (26) numerically. As we 
have infinite energy bands below zero, we introduce a 
cut-off function g(s a ) which smoothly vanishes \e a \ > s c . 
In the following, we take e c = 50eo where 



£o = 



2-kHv 



(27) 



is the typical energy scale for the subband structure. The 
result is actually converging in the limit of large e c . 

The susceptibility of the infinite bulk graphene at zero 
temperature is given by 1,13,20 



The operator of the electric current density is given by 

Jy(r) = -\ {v y S(r - r>) + S(r - r')v y } . (21) 

In the first order perturbation in S~H, the expectation 
value of the current density is written as 



2Re 



{VH-) a pQy{r))fi 



ep 



(22) 



where a and ft represent the unperturbed eigenstates of 
graphene ribbon, and /(e) = 1/[1 + e^^ E ~^] is the Fermi 
distribution function with the chemical potential fx. 

In a zigzag ribbon, the current density always vanishes 
at the edges x = and L x , while it is not generally the 



Xgr(£F) = -g v g s -^r^S(e F ) 



(28) 



where g v — g s — 2 are the valley and spin degeneracies, 
respectively. At finite temperature, it becomes 



Xgr^;? 1 ) 



de 



dm 

de 



= -g v gs 



2 2 

e v 



Xgr(e) 
1 



24ttc 2 fc B Tcosh[^/(2fc B T)] ' 



(29) 



In the graphene ribbon, the characteristic unit of the 
susceptibility can be chosen as 



Xo = g v g s 



e v 



Que 2 £ 



(30) 
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FIG. 3: Band structure (upper panel) and magnetic susceptibility as a function of ef (lower), of (a) zigzag, (b) metallic armchair 
(y — 0) and (c) semiconducting armchair (y — ±1) graphene ribbons. In upper panels, solid (black) and dashed (red) curves 
indicate the band structures at zero and a finite magnetic field, respectively. For the latter, the energy band is calculated with 
the perturbation theory in a magnetic field B. where we take B = Bo for (a), and B — 0.5-Bo in (b) and (c) for illustrative 
purpose, and Bo = (h/e)/L 2 x . 
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FIG. 4: Magnetic susceptibility against chemical potential 
\i for metallic armchair ribbon (y = 0) at several different 
temperatures. 



III. NUMERICAL RESULTS 
A. Magnetic susceptibility 

Lower panels of Fig. 3 show the orbital susceptibility 
X{ £ f) of (a) zigzag, (b) metallic armchair (y = 0) and 
(c) semiconducting armchair (y = ±1) ribbons at zero 
temperature, where the upward direction represents the 
negative (i.e., diamagnetic) susceptibility. The figures 
are to be compared with the band structures in upper 
panels. In every case, the magnitude of x becomes the 
maximum at Sp — 0, and oscillates as a function of Sp 
in accordance with the subband structure. In the posi- 
tive energy region, for example, the curve sharply rises 
when a subband starts to be occupied by electrons, while 
it tends to decrease otherwise. In large |e_f|, the ampli- 
tude of the oscillation slowly attenuates approximately 
in proportional to 1/ \/|£f[. 

The oscillating feature can be understood in terms of 
the band energy shift in an infinitesimal magnetic field. 
In the upper figures of Fig. 3, we plot as broken curves 
the energy band in some small B calculated by the second 
order perturbation. Here the amplitude B is set to some 
finite value for illustrative purpose. Generally the system 
is diamagnetic when the total energy shift caused by B is 
positive, and paramagnetic when negative. In the metal- 
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FIG. 5: Diamagnetic current density j y (x) of different types 
of graphene ribbons with ef = at T = 0. 
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lie armchair ribbon (b), for example, we see that a pair 
of the first subbands (n = 0, s = ±) shift towards zero 
energy, due to the level repulsions from excited subbands 
nearby. All other bands (n > 0) move in the opposite di- 
rection away from zero energy, while the absolute shifts 
are much smaller than that of n = 0. 

When sf is zero, the energy gain of the first va- 
lence band (s = — ,n = 0) exceeds the energy loss of 
all other valence bands (s = —,\ n \ ^ 1)> resulting in 
the total diamagnetism. When ep is shifted to positive 
side, the diamagnetism decreases because the first con- 
duction band (s = +, n = 0) has a negative shift and 
gives paramagnetism. When the second conduction band 
(s = +, \n\ = 1) starts to be filled, the susceptibility sud- 
denly jumps to diamagnetic direction, because the shift 
is positive there and also the density of states diverges at 
the band bottom. The oscillation of other types, (a) and 
(c), can be explained in a similar manner. 

In Fig. 4, the susceptibility at several different temper- 
atures is plotted as a function of the chemical potential 
jjL. We here choose the metallic armchair ribbon {v = 0) 
while the qualitative property is the same in other cases. 
We define the characteristic temperature scale To as 

k B T = e , (31) 

at which the thermal broadening energy is of the order of 
the subband interval energy. We see that the oscillation 
rapidly disappears once T becomes of the order of To, 
leaving only single diamagnetic peak at sp — 0. When 
T >^ To, the curve becomes almost identical with the 
bulk susceptibility, Eq. (29). We also confirmed that the 
integration of \ over sf is identical with the bulk value 
— g v g s e 2 v 2 / (6ttc 2 ) within the numerical accuracy, for all 
the types of ribbons considered here. This fact suggests 
that the finite size effect disappears T >^ To, regardless 
of the edge configuration. 



FIG. 6: Two-dimensional density plot j y (x;EF) of the dia- 
magnetic current density of metallic armchair ribbon (y — 0), 
as a function of position x (horizontal axis) and Fermi energy 
(vertical). 
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FIG. 7: Diamagnetic current density j y (x) of metallic arm- 
chair ribbon [y — 0) with ef = 0, at several different temper- 
atures. The vertical arrows indicate the characteristic length 
scale ftu/(27rfcsT) measured from x = 0. 



B. Diamagnetic current density 

Fig. 5 shows the diamagnetic current density j y (x) in 
different types of graphene ribbons at ep — and T = 0, 
calculated in the first order perturbation of B. The unit 
of current density is taken as c\oB / L x . The current flows 
in opposite directions in the left-hand side and right-hand 
side of the ribbon, to make a magnetization perpendicu- 
lar to the layer. Reflecting the absence of the character- 
istic length scale, the current distribution is not localized 
to the edge but spread in the entire width in a form of 
slowly-varying monotonic function. 
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In zigzag ribbons, the current density actually becomes 
absolute zero at x = and L x , in accordance with the 
constraint argued in the previous section. The current 
sharply drops to zero at the edges, and some oscillatory 
feature remains around the edge due to a finite cut-off en- 
ergy. When we increase the energy cut-off (not shown), 
the curve appears to slowly approach a fixed curve hav- 
ing a discontinuous jump at the edges. In the armchair 
ribbons, j y is not necessarily zero but logarithmically di- 
verges at the both edges. The numerical calculation con- 
verges much more rapidly there, since there is no discon- 
tinuity as in the zigzag case. 

Fig. 6 is the two-dimensional density plot j y (x;£p) of 
the diamagnetic current density of metallic armchair rib- 
bon, as a function of position x (horizontal axis) and 
Fermi energy (vertical). In increasing £p, the current 
distribution begins to oscillate as a function of x, with a 
characteristic wave length of the order of kp — ef/{Kv). 

The temperature dependence of the current density at 
£f = is shown in Fig. 7 for the same metallic arm- 
chair ribbon. When T becomes as large as T , the cur- 
rent distribution is localized at the boundary forming 
the counter edge currents. This is the same temperature 
range where the oscillation of x disappears and the bulk 
limit is achieved. The depth of the current distribution 
is characterized by 



Aedge (T) 



hv 



2irk B T' 



(32) 



which shrinks in increasing temperature. With the band 
velocity of graphene, v w 10 6 m/s, it is estimated as 
Acdgc « [l/T{K)]fjan. 

This behavior is intuitively explained using the plot of 
j y {x;EF) in Fig. 6. The current density at a finite T is 
obtained by integrating j y (x; Ef) in ef with the thermal 
averaging factor —df/ds. The current of the middle part 
of the ribbon vanishes in averaging the oscillating func- 
tion in ef, while the cancellation is not complete only 
near the edges, since j y is always positive and negative 
in left and right ends, respectively. The similar temper- 
ature dependence of the current distribution is found in 
other types of ribbons considered here. This suggests 
that, in any finite pieces of graphene with length scale L, 
the finite-size effect disappears when T > To, and then 
the diamagnetic current circulates only near edge with a 
depth A odge . 



C. Relation to spin paramagnetism 

We neglect the effect of the electron spin through out 
the present analysis. In a zigzag ribbon, particularly, the 
large density of states contributed by the zero-energy flat 
band is expected to give a significant magnitude of Pauli 
paramagnetism and reduce the orbital diamagnetism. 31 

The ratio between two effects can be quantitatively 
estimated as follows. The susceptibility of Pauli param- 



agnetism is given by 

Xpara 



(f) »%D(e), 



(33) 



where g ~ 2 is the ^-factor for graphene electron, /ig = 
eh/(2mc) is the Bohr- magneton with to being the free- 
electron mass, and D(e) is the density of states per area 
given by the zero-energy flat band. Since the number of 
edge states accommodated in a ribbon of the length L is 
~ L/a per spin and per valley, 31 we have 



D{e) 



9v9s 
L 2 



L 



(34) 



which gives a delta-function singularity in Xpara- 

By comparing Xpara with the bulk orbital diamag- 
netism Xdia, Eq. (28), we obtain 

1.0 x^, (35) 



Xpara 



Xdi; 



3tt 



TO 



■v 2 aL 



which is negligible in a wide strip with L ^> a. In a low 
temperature such that ksT <~ £q, however, Xdia cannot 
be regarded as thermally broadened delta-function due 
to the effect of the subband formation, and then Xpara 
overcomes Xdia only at £p = 0. 



IV. CARBON NANOTUBES 

The carbon nanotube is a quasi-one-dimensional 
system similar to graphene ribbon, but different in 
that there are no edges. 42 ' 43 Experimentally, graphene 
nanoribbons with smooth edges can be obtained by un- 
zipping the carbon nanotubes, i.e., lengthwise cutting of 
carbon nanotube side walls. 40 ' 41 Then we may ask which 
of the ribbon and the original nanotube has greater dia- 
magnetism, and how the susceptibility oscillation in ex- 
changes in unzipping. The orbital susceptibility of car- 
bon nanotube was theoretically studied for small Fermi 
energies in the effective mass approximation 23 ' 24 . Here 
we compute full Fermi energy dependence in parallel 
fashion to the analysis for ribbons. 

A carbon nanotube is characterized by a chiral vector, 



L = n a a + n b b, 



(36) 



where the atom at L on a graphene sheet is rolled up 
onto the origin in constructing a tube. The bound- 
ary condition is given by V'a(-Ra) = iPa{Ra + L) and 
V'b(-Rb) = iPb(Rb + L). For the effective mass wave- 
function, it is written as 53 



F K (r + _£,)= cxp (-^„) F K (r), 
F K '(r + L) = cxp (+™A F K '(r). 



(37) 



Here F K = (Fg,Fg) etc., and v is an integer (0,±1) 
defined by 



n a + n b = 3m + v, 



(38) 
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with integer m. 

For K point, the eigenstates are immediately obtained 

as 



Ff nk (r) = 



y^2L X Ly 



i(k n x+6) 



V o J 



£snk v — shV\/k v + /c 2 , 



where 



2tt 

Lr 



0,±1,±2,- 



(39) 



(40) 



and L x — \L\ and L y is length of the carbon nanotube. 
The system is metallic when v = 0, and semiconducting 
when v = ±1. The band structure looks similar to arm- 
chair graphene ribbon's, but the unit of momentum quan- 
tization doubled compared to Eq. (15), leading to wider 
energy spacing between subbands. The energy band for 
K' is obtained by replacing k y by — k y and also v by —v. 

When a uniform magnetic field B is applied perpen- 
dicularly to the nanotube axis, the vector potential can 
be taken as 



, BL X . 2irx 
A (r)=[0,— S m — 



(41) 



We should note that the expression differs from that for 
ribbon, Eq. (18), because the magnetic field perpendicu- 
lar to the tube surface is not a constant, but a sinusoidal 
function in x. Except for that, the magnetic suscepti- 
bility Xtube(eF) is calculated in the same formula, Eq. 
(26). 

The susceptibility of the carbon nanotube is naturally 
related to that of graphene against a spatial varying mag- 
netic field B{q)$mqx with q — 2n/L x = q . When 
we define the g-dependent susceptibility of graphene as 
Xg r (g) = m(q)/B(q), 16 we obtain a relation 



(Xtubc(£F)) V = ^Xgr(qo;S F ), 



(42) 



Here ( ) v represents an average over a phase factor if 
which twists the boundary condition of carbon nanotube 
as ip(r + L) = cxp(2Triip)ip(r) . Physically, the phase fac- 
tor corresponds to threading a magnetic flux of (h/e)<p 
into the nanotube cross section. 23,53 It changes momen- 
tum quantization of Eq. (40) to fc„ = (2ir/L x )(n + ip — 
J//3), and the flux averaging over ip smears the difference 
in v. The factor 1/2 in Eq. (42) enters because the aver- 
age of the squared magnetic field on the nanotube surface 
is B(q Q ) 2 /2. 

At the zero temperature, Xgr(<7; £f) is explicitly evalu- 
ated as 16 



, s gvg s e 2 v l , . 




2k ? 



(43) 



where kp = \ep \/{Kv) is the Fermi wave number and 9(x) 
is defined by 9(x) = 1 (x > 0) and (x < 0). Using Eqs. 
(42) and (43), the flux- averaged susceptibility integrate 
is shown to be 



Xtubc(£F)de F ) = ~ 



-9v9s 



6ttc* 



(44) 



which is exactly half of graphene's, suggesting that the 
susceptibility is effectively smaller in nanotube than in 
ribbon. This is simply because the i?-field component 
penetrating the lattice plane is smaller in the nanotube 
due to its cylindrical shape. 

The susceptibility before taking flux average can be 
calculated in numerics. Fig. 8 (a) shows Xtube(eF) for 
the metallic (y = 0) and the semiconducting (y = ±1) 
nanotubes, together with the flux average. It has an oscil- 
latory behavior similar to the graphene ribbon's, while x 
in \sf\ > frvq/2 completely vanishes after flux average. 16 
In increasing temperature (not shown), the oscillation 
immediately disappears, leaving a single peak regardless 
of v, similar to the graphene ribbon. 

Fig. 8 (a) compares the susceptibility of a carbon nan- 
otube and that of corresponding graphene ribbon un- 
zipped from the same nanotube. Here we chose a zigzag 
ribbon as an example, when the corresponding nanotube 
becomes an armchair nanotube which is always metal- 
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FIG. 8: (a) Magnetic susceptibility of carbon nanotubes. 
Solid (black), dashed (red) and dotted (green) curves are for 
metallic (y = 0), semiconducting iy = ±1), and the flux av- 
erage, respectively, (b) Magnetic susceptibility of the carbon 
nanotube of v = (solid black) and the zigzag ribbon un- 
zipped from the same nanotube (dashed red). 
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lie (v = 0). 53 The oscillation period of the nanotube is 
approximately twice as large as that of the ribbon, re- 
flecting the wider subband spacing. Overall magnitude 
of x is smaller in nanotube roughly by factor 2. The inte- 
grate of susceptibility in ep differs in factor 2 in numerical 
accuracy, in accordance with the above arguments. 



V. RANDOMLY-STACKED MULTILAYER 
GRAPHENE 

The diamagnetism can be made even greater by stack- 
ing graphenes in three dimensions. The recent experi- 
mental technique realizes a novel kind of graphene multi- 
layer in which successive layers are stacked with random 
rotating angles. 45-47 There it is known that the inter- 
layer coupling is significantly weakened and the Dirac 
cone is kept almost intact near zero energy as long as the 
rotating angle is not too small. 48-51 The orbital suscep- 
tibility of such a system is expected to be much stronger 
than graphite in which the delta-function peak of x( £ f) is 
much broadened and shortened by the regular interlayer 
coupling. 18 ' 19 

Here we consider the orbital diamagnetism of a finite- 
sized piece of random-stacked graphene multilayers. In 
calculations, we self-consistently include the effect of the 
counter magnetic field induced by the diamagnetic cur- 
rent itself. This is essential because, as we will show in 
the following, the counter magnetic field of this system 
can be of the same order of the external magnetic field, 
and even nearly perfect screening is possible in low tem- 
peratures. 

For simplicity, we completely neglect the interlayer 
coupling and regard the system as a set of independent 
single layer graphenes. We also assume that each layer 
has the identical shape with a characteristic length scale 
L, and that the system is large enough that the thermal 
broadening energy fc^T is much larger than 2nhv/L. Ac- 
cording to the previous discussions, we then expect that 
the susceptibility of each layer is given by the bulk limit 
Xgi in Eq. (29), and also the depth of the edge current 
•^edge of Eq. (32) can be neglected with respect to the 
system size L. 

Let us consider a situation where a external field B ext is 
applied perpendicularly to graphene plane of the random 
stacked multilayer. The total magnetic field B penetrat- 
ing the system is 



B = B cxt + AS, 



(45) 



where AB is the counter field caused by graphene elec- 
trons. The total field B induces the magnetism in each 
layer, M — x gr i?5, with S being the area of the layer. 
This is related to the diamagnetic edge current I of each 
single layer by 



T CM D 
1 = — = CXgrB. 



(46) 



Since the ring current I exists every interlayer distance 
d, it induces a counter-magnetic field inside the system 
as 



AB 



4nl 
c d' 



(47) 



Solving the set of equations, we find that the dimen- 
sionless volume susceptibility becomes 



X3D 



-1 



AB 



(48) 



At the charge neutral point \i = 0, in particular, we have 

-1 



X3d(a* = 0) = 



1 + fcsT/A' 



where A is a characteristic energy scale defined by 



9v9s 



0.03meV, 



(49) 



(50) 



and d is assumed to be the interlayer spacing of graphite, 
0.334 nm. 

In decreasing the temperature, xw monotonically in- 
creases in the negative direction, and approaches — 1, 
where the prefect magnetic field screening is achieved. 
This reflects the property of the single-layer susceptibil- 
ity, Eq. (29), of which peak value at Dirac point diverges 
in T — > 0. In contrast, X3B of the graphite is of the order 
of 10~ 4 and is not much enhanced in low temperatures, 57 
because x( £ f) is already broadened by the interlayer cou- 
pling energy about 4000 K. 18,19 A three-dimensional bulk 
material composed of random-stacked graphenes, if real- 
ized, would be the strongest diamagnetic material than 
any other known substances except for the superconduc- 
tors. Including the effect of the residual interlayer cou- 
pling between misoriented layers may set the upper limit 
to X3D, while we leave the detailed analysis for a future 
problem. 



VI. CONCLUSION 

We have studied the orbital diamagnetism of graphene 
ribbons using the effective-mass approximation, to figure 
out its dependence on temperature and Fermi energy, 
and also the finite-size effect on the delta-function sin- 
gularity in the bulk limit. In increasing temperature, an 
oscillatory behavior in the orbital susceptibility x( £ f) is 
eventually smeared out, approaching the bulk limit i.e., a 
thermally broadened delta- function centered at ep = 0. 
The electric current responsible for the diamagnetism 
spreads entirely in the sample at T = reflecting the 
absence of the characteristic wavelength, while as T is 
increased, the current density tends to localize near the 
edge with a depth ~ hv / finksT) . 

We also see a carbon nanotube, another form of quasi- 
onc-dimensional carbon, exhibits a similar oscillation in 
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x(sf), but the overall magnitude is reduced by a fac- ACKNOWLEDGMENTS 
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